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1 Introduction 

The inverse problem in EEG and MEG amounts to 
finding a realistic source distribution in the human 
brain for a given set of field observations on the sur¬ 
face of the head. This requires the repeated solution 
of the forward problem, i.e., the simulation of the field 
propagation for a given dipolar source in the brain 
using a volume-conduction model of the head. For 
the most realistic modelling, the different tissues have 
to be segmented and assigned individual conductiv¬ 
ity tensor material parameters. Layers like the skull 
or fibrous tissues like brain white matter are known 
to be anisotropic up to a ratio of 1:10. One of the 
most sensitive parameters is the anisotropic conduc¬ 
tivity of the skull. The identification of the cere¬ 
brospinal fluid (CSF)-skull boundary based on Tl- 
weighted MRI (Tl-MRI) is problematic and Proton- 
Density-Weighting (PD-MRI) is most suitable for this 
task. In this paper we will outline how individual¬ 
ized high resolution finite element (FE) models, ex¬ 
ploiting multimodal MR-imaging protocols, are au¬ 
tomatically constructed. We present an improved 
segmentation of the skull through a combination of 
Tl- and PD-MRI. The structural information about 
the white matter fibre directions are won through an 
MR-diffusion tensor imaging protocol [1]. The use 
of fast techniques to solve the large sparse systems 
of linear equations arising from the 3D FE method 
is necessary in order to have more acceptable so¬ 
lution times with high resolution anisotropic mod¬ 
els and inverse source localization. Preconditioned 
Krylov-subspace-methods are among the most attrac¬ 
tive iterative methods. We will compare incomplete 
threshold-factorization- with multigrid- precondition¬ 
ers, the latter is known to be an optimal method with 
respect to the operation count and memory. Since 
the geometric construction of a grid-hierarchy is diffi¬ 
cult (we only have tensor measurements for the finest 
level), we use a pure algebraic multigrid (AMG) pre¬ 
conditioner. 


2 Methods 

2.1 Segmentation and mesh generation 

The following description summarizes the 
segmentation/registration-results obtained in [2]. 
A voxel-similarity based affine registration without 
pre-segmentation using a cost function based on 
mutual information was used to register the PD-MRI 
onto the corresponding TI-MRI. The maximization 
of the mutual information was carried out using a 
multilevel-downhill-simplex approach exploiting 
Freudenthal-triangulation, which turned out to be the 
most accurate and the fastest method compared to a 
genetic optimization approach or Powells direction 
set method. Initial skull surface segmentations were 
calculated using a fuzzy-C-means classification algo¬ 
rithm which corrects for intensity inhomogeneities in 
the MR-images. Within this procedure, the segmen¬ 
tation of the inner skull surface was obtained from the 
PD-image whereas for the outer skull surface both 
image modalities were exploited. The initial surfaces 
were triangulated and treated as a deformable model 
to obtain the final skull-segmentation result. The 
segmentation of CSF, white and grey matter was 
carried out following [3]. 

Motivated by the results of [4], a voxel-based mesh- 
generator [5] with and without surface-smoothing 
(node-shift) was developed to generate surface- 
smoothed hexahedra FE-meshes which better take 
into account the result of the accurate skull- 
segmentation. 

2.2 Fast solver methods for FE-equations 

A description of the multigrid-method as a precondi¬ 
tioner within a Krylov-subspace method can be found 
in [7]. [8] and [9] focus on incomplete factorization 
preconditioners and their parallelization. 

The condition number of the symmetric positive- 
definite geometry-matrix A of the linear equation sys¬ 
tem A<I> — J, with the given source J and the un- 



Algorithm 1 Setup process for AMG: Setup(A/,l) 
if \u t \ > COARSEGRID then 
Split oji into disjoint sets ooc and 00 f 
Set uiFi = t o c 

Define the interpolation operator P, R — P T 
A 1+1 = RA t P 
Setup(A H i,l+l) 
else 

Perform an LL^-factorization of A[ 
COARSELEVEL = 1 

end if 


known potential <&, is of order k(A) — 0(h~ 2 ), 
where h is the mesh-size of the FE mesh. A large con¬ 
dition number limits the accuracy and performance 
of Krylov-methods so that the goal of a precondi¬ 
tioner, M, is the reduction of n(M~ l A) for the pre¬ 
conditioned equation system M~ l AQ — M~ l J. A 
second requirement is that it is cheap to solve lin¬ 
ear systems Mzj — rj (zy residual for the pre¬ 
conditioned system, ry residual), needed in the j th 
step of the iterative Krylov method. Up to now, 
we used “matrix-pattern” factorization precondition¬ 
ers for the conjugate gradient (CG) method. The 
diagonal scaling or Jakobi-preconditioning chooses 
M — D with D the diagonal entries of A. The 
IC(0)-factorization preconditioner exploits an incom¬ 
plete Cholesky-decomposition M — LL l of A with 
zero fill-in, i.e., L has the same non-zero-pattem as 
A. Such factorizations are “blind” to numerical val¬ 
ues because elements that are dropped depend only 
on the structure of A. In this paper, we tested a 
threshold-technique where elements are dropped ac¬ 
cording to their magnitude, the incomplete LDL 1 fac. 
with threshold value e (ILDLT(e). The performance 


Algorithm 2 V(V F ,^)-cycle: AMG(A/, u h f u 1) 
if l = COARSELEVEL then 

ui — (LL 1 )- 1 fi (forward-back sweep) 

else 

Relax up times on Aiui — fi 
Calculate the defect di — fi — Aim 
Restrict the defect di+i — Rdi 
Set ui +1 = 0 

Apply AMG(i H i,t/| + i,d H i,l+l) 
Prolongate the correction si — Pui+i 
Update the solution ui — ui + si 
Relax up times on Aiui — fi 

end if 


Table 1: Condition numbers of the different models. 


Model 

eij 

largest 

^envalue 

smallest 

condition 

number 

Rea-ns-cube 

14.282 

0.506 * 10“ e 

2.822 * 10 7 

Rea-cube 

13.570 

0.494 * 10~ 6 

2.742 * W 

4Sp-nscu-iso 

8.114 

0.132 * 10 -t) 

6.147 * 10 6 

4Sp-nscu-ani 

8.114 

0.126 * 10“ 5 

6.439 * 10 6 

4Sp-cu-iso 

7.74 

0.13 * 10 -5 

5.85 * 10 6 

4Sp-cu-ani 

7.74 

0.12 * 10 -5 

6.09 * W 

4Sp-tet-iso 

97.28 

0.62 * 10 -5 

1.57 * 10 7 

4Sp-tet-ani 

42.83 

0.56 * 10“ 5 

7.71 * 10 e 


of the ILDLT(e)-CG method will be compared with 
an ILDLT(e) preconditioned symmetric variant of the 
Quasi-Minimal Residual (QMR) algorithm derived 
in [9]. For factorization preconditioners, the system 
Mzj — rj can be solved by a forward-back sweep. 
The multigrid method requires 0(N log r] -1 ) arith¬ 
metical operations to reduce the initial error by the 
factor r] (N: problem size) [6]. The Multigrid- 
preconditioner leads to a condition number of 0(1) 
[7]. The number of iterations is also independent of 
h. For the algebraic multigrid, these result have not 
yet been proved but they are accepted as ‘empirical 
results’. In the following we are concerned with the 
construction of an AMG preconditioner. A^ — A 
can be interpreted as an FE-grid, i.e., the diagonal en¬ 
tries of the matrix Ah are related to grid points in u 
and off-diagonal entries are related to edges in an FE- 
mesh. First we investigate the coarsening process. 
Motivated from an FE mesh, the grid points 00 ^ (or 
equivalently the matrix AC) can be split into two dis¬ 
joint subsets Uh — uoc U( jOf such that there are (al¬ 
most) no direct connections between any two coarse 
grid nodes and the resulting number of coarse grid 
nodes is as large as possible (for more detailed infor¬ 
mation see [10]). Next the prolongation operator P 
from a coarser to a finer grid has to be defined. The 
most simple prolongation and the one which turned 
out to be most efficient in our simulations is the ‘equal 
distribution’ which is given by 

{ 1 i = j <E uc 

f~\ ieuFjeui; 

0 otherwise 

with u> l c = N* fl u>c and N' 1 — {j \ \aij\ > 8 ■ |a*j|}. 





Figure 1: Skull surface segmentation result through 
multimodal MR-imaging presented on underlying 
MR-images. a) T1 -weighted b) registered PD- 
weighted 


Figure 2: a) Axial layer of the node-shifted FE-mesh 
of a realistic 5-tissue model, b) Isopotential-lines 
from —2 pV to +2pV on a coronal layer of the node- 
shifted FE-mesh. 


The restriction operator is defined by R P T . With 
these definitions, the coarse grid operator Am is real¬ 
ized by the classical Galerkin method Ah — RA^P 
and a recursive application of the above steps imme¬ 
diately leads to a multilevel setup (Alg. 1). Finally 
an appropriate smoother is defined by a GauB-Seidel 
method and thereafter a usual multigrid cycle is given 
in Alg. 2. For the m-V(vF,VB) cycle AMG precon¬ 
ditioned CG method, the preconditioning operation 
Mzj — rj is equivalent to m calls of AMG(A,Zj,rjA) 
where the current zj is always used. 

3 Results 

3.1 Segmentation and mesh generation 

The results obtained by the procedure described in 
Section 2.1 are shown in the Fig. 1 and 2a. 

3.2 Solver-comparison 

The isotropic conductivity values of the 5-tissue re¬ 
alistic head model were chosen as follows (in S/m): 
skin 0.33, skull 0.0042, CSF 1.0, grey matter 0.33 and 
white matter 0.14. Isotropic (Rea-cube) and node- 
shifted (Rea-ns-cube) 2mm cube meshes were gener¬ 
ated. A 4-layer sphere model (cond. 0.33, 0.0042, 
1.0 and 0.33) was constructed and meshed with 
2mm isotropic (4Sp-cu-iso) and 2mm node-shifted 
(4Sp-nscu-iso) cube elements. A radial Tangential 
anisotropy of 1:10 was assigned for the ’'skull” layer 
(4Sp-nscu-ani, 4Sp-cu-ani). A geometry-based mesh 
generator was used to generate a tetrahedral mesh 


for the isotropic (4Sp-tet-iso) and anisotropic (4Sp- 
tet-ani) 4-layer-sphere model. After setting up the 
geometry matrices, the condition numbers, defined 
as the quotient of largest over smallest eigenvalue, 
were calculated to obtain an impression about the 
ill-posedness using the algorithms described in [9]. 
The solver-performance tests were carried out on a 
1-processor-machine (CPU: MIPS R10000,180MHz, 
FPU: MIPS R10010, 2MB sec.cache, 32KB instruc¬ 
tion and 32KB data cache). The starting vector for 
all solvers was 0. Since the solver-performance di¬ 
agrams were quite comparable for all tested equation 
systems, only two representative figures are presented 
(Fig. 3 and 4). The selected values of S in the AMG- 
and e in the threshold-case are indicated behind the 
preconditioner, the number of iterations can be found 
behind the Krylov-method. We used the 1-V(1,1)- 
cycle-AMG-preconditioner. The setup-times for the 
preconditioners were neglected since it has to be per¬ 
formed only once for the inverse problem. To give 
an example, the AMG in Fig. 3 used 6 levels with 
sizes 323752, 57929, 8613, 2146, 1126 and 687. The 
setup-time was 80.2 seconds. The setup-time for the 
ILDLT(le-3)-preconditioner was 37.8 seconds. 

4 Discussion 

The subject in Fig. 1 shows an above-average amount 
of CSF. A segmentation of the inner skull surface 
based on only the T1-weighted MRI would be es¬ 
pecially flawed in such cases and could lead to large 
mislocalizations [11]. 
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Figure 3: Solver comparison (L, 2 ~relative residual). 


The results in Table 1 show a minor influence of the 
modelled anisotropy and the node-shift on the con¬ 
dition number. The performance tests with ILDLT- 
and AMG-preconditioning, both taking anisotropy 
into account, showed no remarkable difference in the 
solver times. The calculation times of the Jakobi- 
preconditioned Krylov-solvers were up to a factor 
1.25 longer for the anisotropic models, probably be¬ 
cause of less-clustered eigenvalue spectra. Up to a rel¬ 
ative I^-residual of 10 -3 , the differences between the 
ILDLT- and the AMG-preconditioned Krylov meth¬ 
ods are minimal. From 10“ 5 up to 10 -10 , the AMG- 
preconditioned CG method is the fastest in all tested 
cases and a factor 2 to 3 times faster than the best- 
tuned ILDLT-preconditioned Krylov method. Be¬ 
cause of the loss of about 3 digits for the poten¬ 
tial from the source to the electrodes, the interesting 
residuals begin at about 10 -5 . For very small resid¬ 
uals, the ILDLT-CG shows peaks where the ILDLT- 
QMR has a plateau. This is surely due to the dif¬ 
ferent minimization criteria of both Krylov-methods 
(the CG-method tries to minimize the A-energy-norm 
whereas the QMR method minimizes the i^-norm 
which has been visualized in the figures). In future 
examinations, we will include skull-anisotropy and 
the white-matter tensor measurements, described in 

[1]. Because these tensors change from element to el¬ 
ement, it can be assumed that the eigenvalue spectra 
are more strongly ‘blurred’ than in the tested mod¬ 
els. We conclude that with regard to the inverse prob¬ 
lem, acceptable calculation times can only be reached 
through a parallelization of the presented solvers. 



Figure 4: Solver comparison (L 2 ~relative residual). 
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